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We have carried out calculations of the triple-differential cross section for one-photon double 
ionization of molecular hydrogen for a central photon energy of 75 eV, using a fully ab initio, 
nonperturbative approach to solve the time-dependent Schrodinger equation in prolate spheroidal 
coordinates. The spatial coordinates ^ and t) are discretized in a finite-element discrete-variable 
representation. The wave packet of the laser-driven two-electron system is propagated in time 
through an effective short iterative Lanczos method to simulate the double ionization of the hy- 
drogen molecule. For both symmetric and asymmetric energy sharing, the present results agree 
to a satisfactory level with most earlier predictions for the absolute magnitude and the shape of 
the angular distributions. A notable exception, however, concerns the predictions of the recent 
time-independent calculations based on the exterior complex scaling method in prolate spheroidal 
coordinates [Phys. Rev. A 82, 023423 (2010)]. Extensive tests of the numerical implementation 
were performed, including the effect of truncating the Neumann expansion for the dielectronic in- 
teraction on the description of the initial bound state and the predicted cross sections. We observe 
that the dominant escape mode of the two photoelectrons dramatically depends upon the energy 
sharing. In the parallel geometry, when the ejected electrons are collected along the direction of the 
laser polarization axis, back-to-back escape is the dominant channel for strongly asymmetric energy 
sharing, while it is completely forbidden if the two electrons share the excess energy equally. 

PACS numbers: 33.80.-b, 33.80.Wz, 31.15.A- 



I. INTRODUCTION 

A measurement of the complete breakup of the atomic 
helium target by xuv radiation was achieved over 10 years 
ago Since then, rapid developments in strong xuv 
light sources and momentum imaging techniques have 
made it possible to record all the reaction fragments, 
nuclei and electrons, in double photoionization of the 
simplest two-electron hydrogen/deuterium molecule by 
one-photon absorption [2h5|; and most recently also for 
two-photon absorption Q. For the double ionization 
of H2 by single-photon absorption, only randomly ori- 
ented molecules were investigated in earlier experiments 
(e.g. 01 )• Using "fixed- in-space" molecules, more recent 
experimental efforts include the measurements of energy- 
and angle-resolved differential cross sections by Weber et 
al. for cither equal energy sharing 0, 3 or asymmetric 
energy sharing Q , and by Gisselbrecht et al. @ for equal 
energy sharing at a photon energy of 76 eV. These ex- 
perimental studies were at least partially stimulated by 
the goal of understanding the similarities and differences 
between the hydrogen molecule and its atomic counter- 
part helium. However, all recorded fully differential cross 
sections to date suffer from some experimental uncertain- 
ties regarding the alignment angle of the molecule with 
respect to the polarization vector of the laser and the 
emission angles of the photoelectrons. 

From a theoretical point of view, the hydrogen 
molecule exhibits a significant complexity compared to 
helium and, therefore, provides an enormous challenge 
to a fully ah initio description inherent in a multi- 



center, multi-electron system. The single-center con- 
vergent close-coupling method was used to model the 
double ionization of H2 by Khcifets and Bray Q. 
Later McCurdy, Rescigno, Martin and their collabora- 
tors [lo| - [l2j implemented a formulation based on time- 
independent exterior complex scaling (ECS) in spheri- 
cal coordinates, with the origin of the coordinate system 
placed at the center of the molecule, to treat the double 
photoionization at a photon energy of 75 eV. The radial 
coordinates of the two electrons are measured from the 
center, and the radial parts of the wave function were 
either expanded in B-splines or using a finite-element 
discrete- variable representation (FE-DVR). 

The time-dependent close-coupling (TDCC) 
method [l3j . again in spherical coordinates, was 
also extended to calculate the triple-differential cross 
section (TDCS) for double photoionization of the H2 
molecule. While the agreement between the published 
TDCSs from the ECS [ll!, [13 and TDCC [ll calcu- 
lations is basically acceptable, noticeable discrepancies 
remained for a few particular geometries. In the parallel 
geometry, for instance, where the molecular axis C is 
chosen along the laser polarization vector e, the coplanar 
TDCS predictions from the ECS and TDCC calculations 
differ by up to 40 percent when one of the electrons (we 
will refer to it as the "fixed electron" below) is observed 
along the direction perpendicular to the e-axis. In some 
other cases, there exists a noticeable "wing" structure 
in the published TDCC predictions for equal energy 
sharing. Additional TDCC calculations [IJ] suggest 
that the agreement can be systematically improved, 
albeit the above-mentioned discrepancy still exists at a 
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somewhat reduced level. 

Another independent approach [l5| to this problem is 
the very recent timc-indcpcndent ECS treatment, for- 
mulated - as in the current work - in prolate spheroidal 
coordinates. Quite surprisingly, the results of that calcu- 
lation differed from both the earlier ECS and also the 
TDCC predictions, both obtained in spherical coordi- 
nates. Specifically, the ECS prolate spheroidal calcula- 
tions showed differences from the earlier spherical coor- 
dinate calculation for the TDCSs, at a level of about 
20% depending on the details of the energy sharing. As 
will be demonstrated below, we have gone to consider- 
able lengths in an attempt to resolve these discrepancies. 
However, significant differences between the present re- 
sults and those of the ECS still remain. 

Both the ECS [nl, [13 and the TDCC [13 calcula- 
tions made some attempt to deal with the experimental 
uncertainties in the scattering angles. Given the experi- 
mental uncertainties and the differences in the previous 
calculations, however, it appeared worthwhile to inves- 
tigate the computational effort required to obtain accu- 
rate TDCSs before averaging over any experimental ac- 
ceptance angles. Consequently, the present calculation 
represents an independent implementation of the time- 
dependent FE-DVR approach in prolate spheroidal coor- 
dinates. As in the other approaches mentioned above, the 
internuclear separation (R) was held fixed at its equilib- 
rium distance of 1.4 bohr. The two-center prolate system, 
with the foci located on the nuclei, provides a suitable 
description for the two-center characteristics of the H2 
molecule. The formulation of the Schrodinger equation 
in prolate spheroidal coordinates for diatomic molecules 
is not new. The pioneering work of Bates, Opik, and 
Foots for the ion, which is exactly separable in 
prolate spheroidal coordinates, already revealed the ap- 
pealing features of the prolate system. In particular, the 
electron-nuclear interaction is rendered benign in this co- 
ordinate system. A partial list of recent applications of 
prolate spheroidal coordinates to diatomic molecules can 
be found in [T^^Hl. 

As has been demonstrated in a number of recent pub- 
lications, a grid-based approach provides a very ap- 
propriate description of laser-driven atomic and molec- 
ular physics when combined with an efficient time- 
propagation method such as the short iterative Lanczos 
(SIL) method [13, [l^. In the present work, we employ 
the FE-DVR/SIL approach in prolate spheroidal coordi- 
nates to study the correlated response of a two-electron 
molecule in the double ionization process. 

The remainder of this manuscript is organized as fol- 
lows. After presenting the Hamiltonian of the hydrogen 
molecule in Section |ll] and providing some details about 
the discretization of the system in an FE-DVR basis in 
Sec. IIII[ the solution of the two-center Foisson equation 
is presented in Sec. IIVI This is followed by a description 
of the procedure for extracting the cross sections of in- 
terest in Sec. |Vl The results are presented and discussed 
in Sec. IVIl before we finish with a summary in Sec. IVIII 



II. THE SCHRODINGER EQUATION IN 
PROLATE SPHEROIDAL COORDINATES 

The prolate spheroidal coordinates with the two foci 
separated by a distance R are defined by 
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and the azimuthal angle ip. Here ri and r2 are the dis- 
tances measured from the two nuclei, respectively. These 
coordinates are specified in the ranges ^ S [l,-|-oo), 
77 e [—1,-1-1], and (f G [0,27r]. The volume element is 
dV = (i?/2)3(^2 _ rf)d£^d'qd(p. According to the asymp- 
totic behaviors as ri,r2 — >■ -l-oo, ^ and rj approach 2r/R 
and cos6', respectively, where r and 9 are the standard 
spherical coordinates. Consequently, f is the "quasi- 
radial" coordinate, while 77 is "quasi-angular" . The 
Hamiltonian of a single electron, Hq, {q = 1,2 for the 
two electrons in H2 below), is written as 



Ha 



1) 



1)9^2 



1 



Vq) dlfil 



(2) 



We solve the time-dependent Schrodinger equation 
(TDSE) of the laser-driven H2 molecule (with two elec- 
trons) in the dipole length gauge: 



4*(l,2,t) = 



H1+H2 + - 
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-E{t)-{ri+r2)\^{l,2,t). 

(3) 

Here is the coordinate of the q-th electron measured 
relative to the center of the molecule and ri2 = |ri — r2| 
is the interelectronic distance. Without loss of general- 
ity, we choose the molecular axis along the z axis, and 
the plane formed by the molecular axis and the polar- 
ization vector as the xz plane. Generally, we can de- 
compose the polarization vector into its two components, 
e = cosOnEz + sin^TvCa;, where 9^ is the angle between 
the <; and e axes, and and are the unit vectors along 
the z and x axes, respectively. The dipole interaction is 
therefore given as 



Eit) ■ (n + rs) =Eit) (zi + Z2) cos 
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where the rectangular coordinates x and z and the pro- 
late spheroidal coordinates are related through 
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We expand the wave function for the H2 molecule in the 
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body-frame as 



m\m2 
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Here <^m,m2i'Pi,'P2) = e'(™i'^i+™='^=)/(27r) is the angular 
function, where mi and m2 denote the magnetic quantum 
numbers of the two electrons along the molecular axis. 

Next, nmim2(^i)^i)C2)'72)i) is expanded in a prod- 
uct of normalized "radial" and "angular" {gkiv)} 
DVR bases: 



n 



Note that the basis is not symmetrized with respect to 
the coordinates of the two electrons. Since we always 
begin the calculation with a properly symmetrized initial 
state, however, that symmetry will be preserved in the 
calculation. 

To discretize this partial differential equation, we em- 
ploy the FE-DVR approach for both the ^ and rj vari- 
ables [2l|, [l^ . If wc normalize the DVR bases according 
to 
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respectively, the overall (^, 77) DVR basis is not normal- 
ized with respect to the volume element in the prolate 
coordinate system. This is corrected by defining the two- 
electron basis 
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X (Cl )/j (6 ).9fc )gi iri2)^mim2 (<Pl , </32 ) , 

which satisfies the desired normalization 

(10) 

to expand n,„^m2 (^1, 771, ^2, ^72, i)- Specifically, wc have 

*(i,2,t)= J2 E^"r^(i'2)^i^rw- (11) 

mim2 ijki 

Introducing the normalized (^, rj) DVR basis eliminates 
the complexities of matrix operations related to the over- 
lap matrix at each time, and hence makes the standard 
SIL algorithm directly applicable to study the temporal 
response of the molecule to laser pulses. 



spheroidal coordinates, we have chosen to work directly 
in the FE-DVR basis. This differs from what is usually 
done for atoms in spherical coordinates, where spheri- 
cal harmonics are used for the angular variables and an 
FE-DVR for the radial coordinates. Consequently, the 
boundary conditions in prolate spheroidal coordinates 
require some more discussion. Analyzing the asymp- 
totics reveals that in the region near the boundaries of 
^ = 1 and ?7 = ±1, which correspond to the molecu- 
lar axis, the single-electron wave function behaves like 
(^2 _ i)|m|/2(i _ ,y2)|m|/2^ rpj^jg indicates that the physi- 
cal wave function is finite for |m| =0, whereas it goes to 
zero in the region close to the molecular axis for |77i| ^0. 
More importantly, the behavior of the wave function for 
odd |to| contains a square-root factor, giving a decidedly 
nonpolynomial behavior to the wave function that is im- 
possible to capture in a straightforward fashion using a 
DVR basis. 

The former problem is readily treated by using a 
Gauss-Radau quadrature in the first DVR element for ^, 
where only the right-most point is constrained to lie on 
the boundary between the first and second finite clement. 
The volume element ensures that the integrand is well 
behaved near the end points and makes it unnecessary 
to invoke a separate quadrature for different m values. 
For all the other elements, a Gauss-Lobatto quadrature 
is employed. This allows us to make the FE-DVR basis 
continuous everywhere and to satisfy the |m|-dependcnt 
boundary condition. 

To overcome the nonanalytic behavior of the basis for 
odd m, Bachau and collaborators explicitly factored 
out the (C^ - l)l'"l/2(l - 772)|m|/2 ^^^^ ^cfore the wave 
function was expanded in terms of i?-splines in the dis- 
cretization approach. We have adopted a similar idea 
in our FE-DVR treatment of the two-center problem to 
achieve much faster convergence, as was also done in 
Ref. [2l|. For the case of even \m\, no changes need to be 
made to define the DVR basis, i.e., the normalized basis 
is written as 
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For odd |m|, however, we define the DVR basis as 
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III. THE FE-DVR BASIS 

In our current implementation of the FE-DVR ap- 
proach for the time-dependent wave function in prolate 



Here and are the weight factors related to the DVR 
bases /i(^) and gi{r]), respectively. The goal of using a 
unique set of mesh points, which are Imj-independent, to 
discretize the (^, ry) coordinates has now been achieved 
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in this scheme. The same technique was employed in 
recent calculations of one- and two-photon double ion- 
ization of H2 [iH Hi- III principle, it is also possible 
to introduce the factors (^2 _ ;^)|m|/2(^ _ ^2)|m|/2 ^^^^ 

the DVR bases to circumvent the difficulties related to 
the nonanalytic behavior near the boundary. However, 
this results in an |m|-dcpcndcncc of the DVR bases and 
quadrature points. This, in turn, leads to a number of 
unnecessary complications in the practical implementa- 
tion of the computational methodology. One might argue 
that an jml-dependent discretization procedure could be 
useful for a system in which the magnetic quantum num- 
ber m is conserved. An example is the ion in external 
magnetic fields along the molecular axis [25|. However, 
that is not the situation in the current calculation. 



IV. THE ELECTRON-ELECTRON COULOMB 
INTERACTION IN PROLATE SPHEROIDAL 
COORDINATES 

Similar to the expansion of the electron-electron inter- 
action in terms of spherical coordinates, a counterpart 
exists in prolate spheroidal coordinates [2g] through the 
Neumann expansion 
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where a = R/2. The two nuclei are located at ±R/2 
along the z axis and C>(<) = max(min)(^i, ^2)- Both 

the regular /';'™'(^) and irregular Legendre func- 

tions [27|, which are defined in the region (1,-t-oo), are 
involved in the expansion as the "radial" part, while the 
"angular" part is only related to P;'™'(7y). Note that 
we chose to work in terms of an un-normalized "an- 
gular" basis rather than the usual spherical harmonics. 
The matrix elements of l/ri2 in a traditional basis, for 
example, a i?-spline or Slater-type basis, can be com- 
puted through the well-known Mehler-Ruedenberg trans- 
formation [1^ [2^ . Due to the discontinuous derivative 
along the line of ^1 = ^2 m the Neumann expansion, 
the straightforward computation of the matrix element 
of l/ri2, using the value of this interaction potential at 
the mesh points, is very slowly convergent. We seek a 
more robust representation of the l/ri2 matrix which 
retains both the underlying Gauss quadrature and the 
DVR property of all potentials being exactly diagonal 
with respect to the highly localized DVR basis. 

In the following, we use the simplified notation 

\ijMm1m2) = |/i(6)/j(6)5fc('7l)5f('72)$mim2(¥'ljV52)) 

to denote the basis. Essentially, we need the integral 
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After integrating over (pi and (^2, the matrix element of 
l/ri2 can be written as 
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where the selection rule 
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= mi — jTi'i = 771 2 — 7712 has 
been used. Hence m is uniquely determined for a given 
pair of angular partial waves. Above we introduced the 
reduced (^, rj) integral 



(18) 

= (ijfc£|p|"|(e<)Ql™i(e>)^^r'(^i)^,'"'(^2)H'/fcr). 

Since |?7i| is fixed in the above equation, we omitted it in 

-^ijke ^ (0 ^-iid will do so in the related quantities below 
as well. It is now worthwhile to define the two electron 
densities 12911: 
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After truncating the radial integral to the edge of the 
box, ^max, this yields 
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Here a convention for the volume element was made in 
such a way that, for any function F{^,r])^ we define 
dV^F{^,r]) = d^a^J^^i^"^ - rf)F{^,ri)dr] to simplify the 
notation. Most importantly, the function iY; (^) is defined 
by 



ui{s,)^Qr\o dv^,pA{e,v')priOpriv') (21) 



Instead of evaluating the above integrals directly, we 
solve the differential equation satisfied by Ui{^). As will 
become apparent later, this equation can be shown to be 
the "radial" Poisson equation in the prolate spheroidal 
coordinate system. The differential equations satisfied 
by the Legendre functions P;'™'(^) and suggests 



that we introduce the operator 
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for given quantum numbers I and m. This is the one- 
dimensional Laplacian operator in the ^ coordinate. 
After some algebra, we obtain 

=^^^^'dy*p^(t,r)pl'"l(i)pl"l(r) (23) 



and 
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Here the Wronskian of the Legendre functions -P/'™'(0 



and (3i™'(0 is given by 

«'(p.""'«).or'K))^{^™- (25) 

Consequently, we obtain 

V|WK6 = e(0- (26) 

This is the second-order inhomogeneous Poisson equation 
satisfied by Ui {£_) with the "source" term given by 
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After carrying out the integral over rj' via Gauss quadra- 
ture, we recast the source term as 

(i~\m\y- (28) 

As one might expect, the inhomogeneous equation re- 
duces to the homogeneous one \i k ^ k' . The solution 
to the Poisson equations and (|28p can be uniquely 
determined by enforcing the boundary conditions 



at ^ = 1 and 



(29) 



Ul{U..) = Ql"'(^max) / dV^,pA{^.^)P]^\OP]^\^) 

(30) 



at C = Cmax, respectively. 

There are two important points to realize, namely: 
First, on the right-hand boundary, the function W/(^) as- 
sumes a nonzero value, which is given by Eq. pop , for all 
possible |m| values. Its asymptotic behavior relies on the 

function (5i"''(0 at large ^, which behaves like l/^fj^ax- 
This indicates that it is nonzero generally, although it 
could be small at the large ^max values used in practical 
calculations. Second, on the left-hand boundary, the sit- 
uation depends on the value of |m|. Z^z(^) takes a nonzero 
value if |m| = 0, while it becomes zero if \m\ ^ 0. 

Following the philosophy employed to handle the 
spherical case [2^, we first seek a solution, Z^°(C), to the 
Poisson equations (P5| - ([5n|) that satisfies the zero-value 
boundary condition at ^max by using exactly the same 
^ mesh points as those for the wave functions. In other 
words, we have Vp"iO = QiO with ^"(1) = Ui{l) and 
^°(Cmax) = 0. After substituting the DVR expansion 

— c^/a'(0 of the solution into the differential 
equation, we obtain a system of linear equations for the 
unknown coefficients {c^}: 



(31) 



><{^1^4)P]^\m). 

The matrix T is defined by its elements 
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Therefore, the coefficient can formally be written as 



[tI' 



^(-l)l"l|[^«^^..'4.'(C^-'7.^)/^^'('7.), 

(33) 

where [pl™!] ^ denotes the inverse of the matrix pl"'l. 
In this case U^{0 fulfills the left-hand boundary condi- 
tion Ul'{l) = 0, and so docs Ui{l). Recall, however, that 
its right-hand boundary condition differs from those of 
W;('Cmax)- This suggests that the final answer to the func- 
tion Ui{^) can be constructed as VliiS,) = ^°(C) + Pi{Oj 
I.e., we add the difference function P,(^) to ^"(0- The 
function P;(^) is also a solution to the homogeneous 
Poisson equation, subject to the boundary condition 
P/(l) = and P/(Cmax) = Z^;(Cmax)- After writing it as a 
linear combination of P;'™'(^) and and imposing 

the boundary conditions, P;(C) takes the form 

F,iO =4.<5,,,a3(ef-r/D^'^'(C0^i"'('7O 

„Ql"'(gmax) (34) 
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We finally arrive at 

m) = 
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|m| = 0. Substituting Eq. (|35]) into Eq. allows us to 
obtain the kernel integral, 



j-i'fk'i' f,^ _ 
-ijke 
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At this point, the DVR version of the solution Ui{S_) is 

ready for all possible values of |m|. cither \m\ 7^ or The matrix clement of l/ri2 can finally be written as 



{ijMmim2\ — \i'i'k'l'm\m'2) ^5u>6jj'5kk'Su'a^i^i ~ vDi^, 
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where we truncated the / summation in the Neumann 
expansion to ^max- The above equation can be converted 
to the normalized (^,77) bases with the help of Eq. ([9|). 
This results in a diagonal representation of the matrix el- 
ements of the electron-electron Coulomb interaction and 
thus considerably simplifies the FE-DVR discretization 
procedure. The above treatment of the l/ri2 matrix was 
successfully applied to the two-photon double ionization 
of H2 [iSl- The implementation of this representation 
will be illustrated below. 



V. TIME EVOLUTION AND EXTRACTION OF 
CROSS SECTIONS 

The time-dependent laser-driven electronic wave 
packet in the hydrogen molecule is obtained by solving 
the TDSE on the {£,,r]) grid. Launched from the pre- 
viously determined ground state, the time evolution of 
the system is achieved by using our recently developed 
SIL method [13, HH. The ground state is determined 
by relaxing the system in imaginary time from an ini- 
tial guess of the wave function on the grid points. At 
each time step we only need to generate the values of the 
discretized wave function on the selected grid points. If 
desired, the information at arbitrary points within the 
spatial box can be obtained from the interpolation pro- 
cedure in terms of the DVR bases. 

A few remarks seem appropriate regarding the efficient 
implementation of the SIL algorithm. The highest en- 
ergy, i?inax, which essentially depends on the smallest 
separation between the (^, 77) mesh points and also on the 



maximum values of |mi| and \m2\, determines the largest 
time step At for the propagation in real time. Typically, 
Euia.K is about 6, 000 atomic units (a.u.) in our calcu- 
lations. Although the chances of electrons populating 
states with such high energies are practically negligible 
for short time scales of the laser-molecule interaction, we 
generally require At < 27r/£'max in order to resolve the 
most rapid oscillations in the time evolution. This means 
that at least a few time steps are needed during one pe- 
riod of 27r/i?inax- Wc refer readers to Refs. [30l - l32l | for 
further details and discussions behind the SIL method. 

In order to ensure that the double-ionization wave 
packet is sufficiently far away from the nuclei, and also 
that the two photoelectrons are well separated, we al- 
low the system to evolve for a few more cycles in the 
field- free Hamiltonian, i.e., after the laser pulse has died 
off. This is the wave packet we use to extract the phys- 
ical information. The ionization probabilities and the 
corresponding cross sections arc extracted by projecting 
the time-dependent wave packet onto uncorrelated two- 
electron continuum states satisfying the standard incom- 
ing boundary conditions. The latter states of H2 are con- 
structed from the one-electron continuum state of the 
ion described in the following subsection. 



A. Continuum states of Hj 

The field-free wave function (f>(^, r/, ip) of the one- 
electron molecular ion is completely separable in pro- 
late spheroidal coordinates. For a given, and conserved. 
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magnetic quantum number m, the wave function takes 
the form $(^,77,(/3) = Tm{0'^ni{il)^ni{'fi), where the az- 
imuthal dependence, is given by ^m{'p) = e™"''/-\/27r. 
The "radial" part Tmq {£,) and the "angular" part Emq {r}) 
of the wave function satisfy the equations 



d_ 



ie-i) 



d 



2i?C + c2^2- A 



and 



— {I- if) — 
drj di] 



(1 ~ rf) 



cW + A. 



771 q 



TraqiO = 

(38) 



(39) 

respectively. Here c = kR/2 for the continuum state 
whose momentum vector has the magnitude k. In addi- 
tion, we need to introduce another quantum number q, 
which denotes the number of nodes of ^miv) the re- 
gion 77 S [—1,-1-1], to label the states, and finally the 
separation constant A^q- 

When the angular function 'E.m{'i])^rni'^) is discretized 
in terms of the relevant DVR bases, a few "spurious" so- 
lutions might be encountered. This is caused by the resid- 
ual errors associated with the Gauss quadratures. Conse- 
quently, we expand the angular function, or "spheroidal 
harmonics" function ytmiViV') = ^mqiv)^mqiyy) with 
£ = Iml+q instead in terms of spherical harmonics. These 
functions arc normalized according to 



d?7 / dipy^^{r],(p)yi,m'ir]^'P) = Smm'Su'- (40) 
1 Jo 



After obtaining the separation constant Amq by solv- 
ing Eq. ((39l) . the "radial" function Tjnq{£,) is again ex- 
panded in terms of the DVR bases. The last DVR point 
at 1^ = ^max needs to be kept for the continuum state. 
Asymptotically, the radial function behaves like 



T, 



1 



-ln(2cO-y 



(41) 

as ^ — > 4-00. Here /S.„iq{k) is the two-center Coulomb 
phase shift. The normalization factor on either the en- 
ergy or the momentum scale and the Coulomb phase shift 
can be determined by matching the numerical solution 
of Tmqi^) according to its asymptotic behavior given in 

Eq. dni)- 

The plane wave in prolate spheroidal coordinates can 
be written as [s^ 

£m 

where i?^^(C) ^ l/(cC)sin[c^ - £7r/2] in the asymp- 
totic region. Note that r]k and rjr are related to the 
directions of k and r in spherical coordinates through 
Vk,r = COS Ok^r- The partial- wave expansion of the plane 
wave e*'^''' reminds us that the two-center Coulomb wave 



satisfying the incoming boundary condition can be ex- 
panded as 



, +00 

X yLik)yi^{Vr,Vr)Tj^^{0- (43) 
This function is normalized in momentum space accord- 
ing to {^'j^^l^k'^) = S{k — k'), provided the asymptotic 
solution in Eq. (|4T|) is satisfied. 

Uncorrelated two-electron continuum states with to- 
tal spin angular momentum S* (S* = in our case) can 
generally be constructed as follows: 

*U.(n,r2) = 



(44) 



1 

71 



s<f,i-)i 



With the help of Eq. P5)) . its partial- wave representation 
can be written as 



*U.(n,r2)= - 



2\3 1 1 



i?/ \/2kik2 



E 



+l2 



(45) 



ijki 



C|^feT'''"'(fcl'fc2) + (-l)^(fcl O fc2) 



Here we introduced 



(46) 



p(fel) (f ^rj.yK2) ( ^■^['^2) / 

il\mi\\'^-^)-^i2\m2\^'^3 '^il\mi\yi'<:)^i2\m2\yl^'^ 

by representing the radial and angular parts on the (^, ry) 
grid points: 



T, 



n(fc2 



,(fc2) 



?^ri(o-E/^(^)^/™(^')' 



^Liv) = E-9A'('?)"£m('?A') 



Here the exchange symmetry 



(47) 

(48) 
(49) 



is satisfied. 



B. Extraction of double-ionization cross sections 

It has been demonstrated 0, HO, [13 that using uncor- 
related two-electron continuum states is a good approxi- 
mation in a time-dependent propagation approach, pro- 
vided the two ejected electrons arc well separated from 
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each other. The probabihty amphtudc of double ioniza- 
tion is then given by 

(50) 



kik2 



X >'^imi(fcl)3^£im2(fe2)5'fimif2m2(^l) ^z), 



where 



Seimie2m2{kl,k2) 
ijki 



(51) 



f2'™2* 



Here the exchange symmetry 

Si2m2limiik2,ki) = (-1)^ deimit2m2{kl, k2) (52) 

is satisfied. We also see that the probability ampli- 
tude formulated in prolate spheroidal coordinates takes 
a similar form as for the atomic case in spherical coordi- 
nates. However, a subtle difference from the atomic case 
is worth pointing out. Strictly speaking, the spheroidal 
harmonics involved in the probability amplitude gener- 
ally depend on the magnitude of the momenta ki and 
k2, in addition to their directions. 

The energy sharing of the two photoelectrons can 
be specified by introducing the hyperangle a — 
tan~^(fc2/A:i). This describes the double-ionization re- 
action with kinetic energies Ei = i?oxc cos^a and E2 ~ 
Eexc sin^a for the two ionized electrons, respectively. 
Here i?cxc is the available excess energy above the double- 
ionization threshold. In the present work, specifically, 
i?exc = 23.6 eV for absorption of a 75-eV photon. 

For double ionization by one-photon absorption, the 
triple differential cross section with respect to one of the 
kinetic energies and the two solid angles fci and fc2 can 
be extracted usin g th e same formalism as in the corre- 
sponding He case [SOtls^]. i.e.. 



dV 



1 



1 



dEidkidk2 



kik2 cos^a Iq t'^-* 



dfc; 



dk' 



X k[S{k2 — k[ tana) 



E 



1+^2 ji^U 



X yiirni (ki , fcl)3^£2m2 (^2 J ^2)'^ Iimil2m2 (^1 J ^2) 



(53) 

Here w and /q are the central photon energy and the 
peak intensity of the laser pulse, respectively, while T^^ 
denotes the effective interaction time between the tem- 
poral electric laser field and the electrons in the one- 
photon absorption process. For a laser pulse of time du- 
ration T with a sine-squared envelope for the field ampli- 



tude, T^l' = (3/8)r. Note that T^^' corresponds to the 



(1) 



special case of the generalized A''-photon effective inter- 
action time Tgfj''' [35l | for a one-photon reaction. Gener- 
alized cross sections for two-photon double ionization of 
the hydrogen molecule were extracted in the same for- 
malism [2J]. 

In the TDCC treatment [l^, different strategies were 
employed to describe the linear one-photon and the non- 
linear two-photon double ionization processes of atoms 
and molecules. For the one-photon case, the cross sec- 
tions were obtained through the time derivative of the 
double-ionization probability, dP^^{t) / dt. The laser 
field does not need to be turned off in this case. On the 
other hand, a true laser pulse was used for the two-photon 
case and an effective time, defined as the time integral 
under a flat-top pulse with a smooth turn-on and turn- 
off, was introduced. In the present work, we employed a 
unified formulation through an effective interaction time 
for both one-photon and multi-photon ionization in laser 
pulses. 

For the one-photon double ionization initialized from 
the lowest X ^Sg state, the two ejected electrons can only 
populate the final and ^n„ continuum states, with 
the specifics depending on the relative orientation of the 
molecular axis and the laser polarization vector. Con- 
sequently only partial waves with ungerade parity [i.e., 
(_l)«i+<?2 = _i] need to be included in Eq. 



VI. RESULTS 
A. Preparation of the initial electronic X^Y^g state 

For the nonsequential double-ionization process in- 
duced by one- or two-photon absorption, electronic cor- 
relation plays a dominant role, as the two photoelectrons 
must share the available excess energy i?oxc- Double ion- 
ization by a single photon would not occur at all if the 
two-electron atom or molecule were approximated by an 
independent-electron model. Therefore, the quality of 
the description of electron-electron correlation in a laser- 
driven system is crucially important for accurate results 
to be obtained. The Coulomb interaction between the 
two electrons has to be described in a consistent man- 
ner for both the initial bound state and the time-evolved 
wave packet. Before we go any further, it is worth dis- 
cussing how we prepare the initial X^T^g state at the 
equilibrium distance of i? = 1.4 bohr. 

As seen from Eq. (|37p . the magnetic quantum num- 
ber in the Neumann expansion of the matrix element of 
l/ri2 is uniquely determined by the angular bases. How- 
ever, this is not the case for the index /, if we choose to 
discretize the coordinate 77, rather then expanding that 
part of the wave function into spherical harmonics. In 
practice, the summation over I must be truncated at a 
finite value of Zmax- In principle, the higher-order ex- 
pansion terms always guarantee well-converged results. 
However, as mentioned earlier, we approximate the rel- 
evant 77-integrals by using Gauss-Legendrc quadrature. 
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This is the price wc have to pay for making the dielee- 
tronic Coufomb potential diagonal in the DVR bases. As 
a consequence, we need to determine how the approxi- 
mation introduced in the 77-integrals for the two-electron 
integrals affects the results for the cross sections of inter- 
est. 

To answer this question, we first investigate the de- 
pendence of the energy obtained for the initial X ^Y,g 
state on the value of ^max used in the Neumann expan- 
sion. Figure [T] shows the variation of the initial-state 
electronic energy of the hydrogen molecule with respect 
to ^max! obtained with a ^ setup of ten elements in the 
region of 1 < ^ ^ 15.82, five in the region 1 < ^ ^ 5 
and another five in 5 ^ ^ ^ 15.82. Each element, in 
turn, contains five DVR points to further discretize the 
configuration space. Furthermore, we employ 9th-order 
DVR points for rj. For a given number of 77 mesh points 
(n^) and |m|i„ax = Iwilmax = |'7i2|max, wc obscrvc that 
the resulting energy typically exhibits a plateau-like be- 
havior with increasing Zmax- For given n^, when ^max 
is relatively small, the 77-intcgral can be computed very 
accurately by using Gauss quadrature. However, when 
^max is too large, the numerical errors introduced from 
the Gauss quadrature cause the energy value to fluctu- 
ate. This occurs when /max approaches 2n^ and is shown 
by the grey stripes in Fig. [TJ In this region of /max an 
unphysically low energy can be produced. Beyond that 
point, the calculated energy increases to the next plateau. 

Ultimately, this is not too surprising, since any Gauss 
quadrature is only reasonably accurate up to a limited 
polynomial order of the integrand. Consequently, if we 
want to keep more terms in the Neumann expansion, 
we have to increase n,, correspondingly. This finding is 
further substantiated by the dependence of the energy 
found for = 11 and 13. The plateaus are indeed ex- 
tended to the correspondingly larger values of 2n,,. Most 
importantly, the amplitude of the energy fluctuation is 
systematically reduced with increasing n^. The error in 
the energy is lowered from 2.13 x 10~^ to 1.45 x 10"'^ 
and finally 1.05 x 10^'^ a.u., when increases from 9 
to 11 and then 13 for |m|max = 4. We obtained the en- 
ergy at i? = 1.4 bohr as -1.8887324 a.u. for /max = 10, 
riri = 9, and |m|max = 4, resulting in a double-ionization 
potential of 51.394 cV. Keeping the other parameters un- 
changed, wc obtained an energy of —1.8887128 a.u. for 
= 11. The benchmark energy in the literature is 
— 1.888761428 a.u. at the same R [3g|, after we take out 
the nucleus-nucleus interaction of 1/1.4 a.u. 

To summarize: Unlike for other expansion parameters, 
it is important to be consistent in the size of the angular 
quadrature and the largest /max employed in the Neu- 
mann expansion in practical calculations, if we discretize 
the coordinate rj. However, this provides a way to ex- 
amine a potential sensitivity of the physical observables 
of interest (here the differential cross sections) to the 
ground-state wave functions generated by varying /max 
and other parameters. This will be further discussed be- 
low. 




10 20 30 40 

^max 

FIG. 1. (Color online) Energy of the lowest electronic X ^Eg 
state at f? = 1.4 bohr as a function of the /max value used in 
the Neumann expansion of l/ri2. The number of 77 points, 
riyj, and the largest magnetic quantum number, Imlmax, are 
labeled as (n,,, |m|max). The open symbols correspond to 
iTilmax = 1, while the filled symbols are for |7Ti|max ~ 4. The 
benchmark energy {Egs) from Ref. [S^ is shown as well. 

B. Convergence of the TDCS 

Before we present our results for the cross sections, let 
us take a closer look at the survival probability 

Psurv = |(*gs|*(t))P (54) 

of the aligned H2 molecule in its ground state ^^gs . This 
is shown in Fig.[2j For homonuclear molecules, the inde- 
pendent alignment angle 0jv between the molecular axis 
and the polarization vector can be confined to the region 
from 0° to 90°. In the xuv regime, wc observe that the 
hydrogen molecule shows a larger probability of being 
ionized or excited (i.e., a lower probability of staying in 
the initial state) at the end of the pulse in an aligned 
geometry. This indicates that the perpendicular compo- 
nent of the temporal electric field exerts more influence 
on the ionization process due to the larger dipolc mo- 
mentum. Interestingly, at the earlier stages of the time 
evolution (e.g., t < 9 a.u.), when the ionized wave packet 
is driven back by the change in direction of the electric 
field, the tilted molecule has a larger probability of stay- 
ing in its ground state. This happens near the various 
minima in Psurv However, once the electric field has be- 
come sufficiently strong {t > 9 a.u.), the wave packet is 
driven out and spread into a larger space. This leads to 
lower minima in Pgurv for the tilted molecule. 

When the wave packet is driven back to the nuclear 
region and therefore has a chance to recombine with the 
HJ ion, a maximum in Pgurv appears. Not surprisingly, 
the parallel geometry always has the largest probability 
for this to happen. Although the wave packet can also 
be scattered for the untilted molecule in the plane per- 
pendicular to the molecular axis, the probability is un- 
doubtedly larger if the laser electric field is perpendicular 
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to the molecular axis. A similar behavior of in xuv 
pulses was observed in Ref. (STj . 
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FIG. 2. (Color online) Survival probability of the hydrogen 
molecule subjected to a sine-squared laser pulse with a peak 
intensity of 10^^ W/cm^. The laser pulse lasts for 10 optical 
cycles and the system is followed for a period of another 2 cy- 
cles of field-free propagation. The central photon energy of 
the laser pulse is 75 eV. 

For most calculations performed in this study, we ex- 
pose the hydrogen molecule to a laser pulse with a peak 
intensity of 10^^ W/cm^. Looking at Fig. [2] we see that 
the depletion of the initial ground state can be safely 
neglected for our typical interaction times. Even for 
9n = 90°, Psurv — 0.99677 remains very close to unity. 
The negligible depletion of the ground state suggests 
that the concept of cross sections is valid and applica- 
ble. On the other hand, it also presents a numerical 
challenge to predict the cross sections accurately from 
a time-dependent treatment, due to the generally small 
ionization probability. 

At first glance, a peak intensity of 10^^ W/cm^ might 
seem very intense for most atomic and molecular targets. 
Here, however, we consider an xuv rather than an IR 
pulse. For an xuv pulse with central photon energy of 
75 eV, such laser fields definitely fall into the "weak-field" 
regime. The ponderomotive energy in the xuv regime is 
much smaller than the photon energy of interest. 

In this work, we are mainly interested in the triple- 



TABLE I. The discretization and expansion parameters of the 
II2 wave function in prolate spheroidal coordinates. Here 
stands for the border between the inner and outer regions 
in the ^ coordinate, while ^max is the size of the ^ box. In 
addition, denotes the number of ^ mesh points in each 
element. The numbers of ^ elements in the inner and outer 
region are ninn and Tiout, respectively. These ^ parameters 
produce the total number of ^ mesh points N^. The ^ grid I 
and 5 grid II are used to examine the convergence of our 
results. 
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^inn 




^max 




N, 




grid I 


5 


5 


67 


150 


5 


288 




grid II 


9 


1 


11 


97 


14 


156 



differential cross section, since it reveals the fine details 
of possible energy sharings and preferred directions of the 
ejected electrons in the double-ionization process. Given 
the discrepancies between results from various calcula- 
tions found in the literature, we carried out comprehen- 
sive convergence tests for our predictions of the TDCSs. 
These tests are essentially divided into two groups. The 
first group concerns the laser parameters, while the sec- 
ond one deals with the discretization and expansion pa- 
rameters. An example of two different parameter sets for 
the ^ grid is given in Table |T] and will be further discussed 
below. 

In order to obtain a good handle on the sensitivity of 
the results to the various parameters and the resulting 
level of "convergence" , we try to only vary a single pa- 
rameter while keeping all others fixed if possible. For 
the dependence on the laser parameters, we use the ^ 
grid I combined with (n^, |m|niax, 'max) = (9,4, 10). For 
the tests regarding the discretizations and expansions, 
the peak intensity of laser was fixed at 10^^ W/cm^ and 
a time scale of "10 + 2" optical cycles (o.c.) was used. 
Here "10 -I- 2" refers to a 10-cycle laser pulse with a sine- 
squared envelope for the field amplitude, followed by a 
2-cycle field-free propagation. 

Figures [3] and |4] show the convergence pattern of our 
TDCS results for asymmetric energy sharing in the par- 
allel geometry {9^ = 0°). The energy sharing between 
electron 1 (observed at the fixed angle 9i) and electron 2 
(observed at the variable angle 02 ) is 20% : 80%. Only 
the electron that takes away 20% of the excess energy is 
recorded at fixed positions cither parallel or perpendicu- 
lar to the polarization axis. 

Since the laser pulse is explicitly involved in our time- 
dependent treatment, we first have to be sure that the 
extracted cross sections are essentially independent of the 
laser intensity and the time scales. Only then are the cal- 
culations of cross sections meaningful. This also allows 
us to compare the physical information extracted from 
our time-dependent scenario to that obtained through 
conventional time-independent treatments, which are ef- 
fectively equivalent to the weak-field approximation and 
"infinitely" long interaction times. Rather than comput- 
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FIG. 3. (Color online) Convergence of the coplanar TDCS 
results for the hydrogen molecule for asymmetric energy shar- 
ing with respect to the laser peak intensity and the time scale. 
The central photon energy is 75 eV. The slow reference elec- 
tron, observed at the fixed angle 6\, takes away 20% of the 
available excess energy [Ei = 4.7 eV), while the other electron 
takes 80% of iSexc {E2 = 18.9 eV). The peak laser intensity 
in panels {b)-{d) is 10^^ W/cm^. The two columns show the 
corresponding convergence of the TDCS for Oi = 0° (left) and 
01 = 90° (right), respectively. 1 barn (b) = 10"^* cm^ 



ing the cross sections, it would be more appropriate to 
consider ionization rates if the cross sections were found 
to be sensitive to the laser parameters. 

In Fig. [21 we display the dependences of our TDCS re- 
sults upon the laser parameters. Note that the TDCSs 
extracted from /o = lO^^ W/cm^ and 10" W/cm^ at 
fixed time evolution of "10 -I- 2" cycles are nearly identi- 
cal and agree with each other to better than the thick- 
ness of the line. When we turn to the dependence of 
time scales at a fixed intensity of 10^^ W/cm^, we use 
the same pulse, but allow the system to freely evolve for 
a few additional cycles to extract the TDCS. This corre- 
sponds to the time scale of "10 -I- 4" o.c. Also, we may 
increase the laser-molecule interaction time, but extract 
the TDCSs at the same cycles of field-free time evolu- 
tion after the pulse died off. This gives the scenario of 
"12 -I- 2" o.c. Since the total time durations are the same 
(14 o.c), they allow us to examine the extracted TDCSs 
from different perspectives. The increased interaction 
time yields a reduced bandwidth of the photon energy, 
while the longer field-free propagation ensures that the 
double-ionization wave packet is further away from the 
nuclear region [38|. The calculated TDCSs indeed show 
a slight, though in our opinion acceptable, sensitivity to 
the time scales. Not surprisingly, the sensitivity is most 
visible for the smaller cross sections, when the two ejected 
electrons travel nearly parallel along the same direction 
(c.f. Fig. He)). 

Having confidence in using the current sets of laser 
parameters, we now turn our attention to the scheme of 
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FIG. 4. (Color online) Same as Fig.O but for the convergence 
of the TDCS results with respect to the discretization and 
expansion parameters. See text for the details. 



spatial discretization {n^, ^max, 't-??) and the convergence 
of the expansion (/,„ax, |™|max)- The results are displayed 
in Fig. m 

For the discretization parameters, we obtain well- 
converged TDCSs by increasing from 5 to 7, n,, from 9 
to 11, and extending the spatial box of ^max from 100 to 
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FIG. 5. (Color online) Coplanar TDCS of the aligned hydrogen molecule for equal energy sharing [Ei = E2 — 11.8 eV). The 
central photon energy is 75 eV. One electron is detected at the fixed direction of 9i — 0° with respect to the laser polarization 
axis. Also shown are the one-center spherical ECS results [ll|], the two-center prolate spheroidal results [l^, and one-center 
spherical TDCC results 



150. Most importantly, however, we consider two sets of 
^ mesh points: ^ grid I and grid II (see Table |T|. The 
principal motivation was to sec whether or not we can 
reproduce the much lower TDCS values (by about 20% 
compared to the one-center spherical results) that were 
recently obtained in an ECS calculation in two-center el- 
liptical coordinates by Tao et al. [T5| . 

We emphasize that these two grids in the "radial" ^ co- 
ordinate are completely different regarding both the dis- 
tribution of the elements and the number of grid points 
per element. In the ^ grid I, we divide the ^ space into 
two parts, an inner and an outer region with a border at 

= 5. We place a narrow span of elements in the inner 
region, and then wider elements in the outer region. In 
contrast to that, ^ grid II does not distinguish between 
inner and outer regions, i.e., the elements uniformly span 
the region from 1 to ^max- The mesh setup in ^ grid II is 
the same as that used in Ref . [l^ , except for the complex 
rotation. The ^ grid I has a much denser distribution of 
mesh points than ^ grid II. Nevertheless, the extracted 
TDCSs from both sets of ^ grids are in excellent agree- 
ment with each other, even for the smallest cross sections. 



This strongly suggests that the results are well converged 
at least with regard to the ^ grid. Both ^ sets are good 
enough to capture the physics of interest. Differences at 
the 20% level are unlikely to be caused by using different 
sets of ^ meshes. 

Finally, we discuss the convergence of our results with 
respect to the expansion parameters, | 

^|max and /max- 

As expected for a onc-photon process, |TO|max — 4 pro- 
duced well-converged results. 

Recall the discussion above regarding the ground state, 
especially how the truncated Neumann expansion of 
l/ri2 in our present FE-DVR implementation affects the 
initial-state energy and therefore the quality of the wave 
function. For consistency, we use the same Zmax iii the 
real-time propagation and in the ground-state wave func- 
tion. As seen from Figs.|4|i) andlU^j), our truncated Neu- 
mann expansion has little effect on the calculated TDCS 
values. Well-converged results can be obtained even with 
an inappropriately large value of /max = 20, which yields 
a slightly higher energy of the ground state (c.f. Fig. [1]). 

Overall, our detailed convergence tests only reveal a 
very weak sensitivity of the TDCS results to both the 
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FIG. 6. (Color online) Same as Fig. [5] except that the fixed electron is detected at the angle 6i = 90° with respect to the 
laser polarization axis. Since there was a plotting error in Fig. 3 of Tao et al. [31, we are comparing here with the proper 
numbers [39| from that calculation. 



time scales and the values of Zmax- Well-converged TDCS 
results can be obtained by using either ^ grid I or ^ grid II 
combined with {rir,, |m|max, ^max) = (9, 4, 10). In the pro- 
duction calculations for the TDCSs shown in the next 
subsection, we used the ^ grid I to discretize the two- 
electron wave packet and a "10 + 2" sine-squared laser 
pulse with a peak intensity of 10^^ W/cm^. 



C. TDCSs for the aligned H2 molecule 

Figures [SJ [6l and [7] display the coplanar TDCSs of 
the aligned hydrogen molecule at equal and asymmet- 
ric {El : E2 = 20% : 80%) energy sharing. The two 
electrons are detected in the same (coplanar) plane de- 
fined by the ^ and e axes. The angles 9i, 62, and Ojs; 
are all measured with respect to the laser linear polariza- 
tion axis. We compare our TDCS predictions with those 
obtained in the time-independent one-center spherical 
ECS calculation [ll|, the time- independent two-center 
spheroidal ECS model [l^, and the time-dependent one- 
center spherical TDCC approach The TDCC num- 
bers were recently recalculated with a bigger box size 



and differ, in some cases substantially, from those pub- 
lished originally [l^. Except for the recent two-center 
prolate spheroidal ECS results of Tao et al. [HI, H^, the 
agreement between the other three sets of results is very 
satisfactory. Once again, the largest relative differences 
occur when the cross sections are small (see Figs, ^d) 
andinid)). 

Using spheroidal coordinates as well, as an illustrative 
example of their two-center ECS approach, Serov and 
Joulakian [i^ recently presented the TDCS at the same 
photon energy, but only for a single geometry of 9n = 20° 
and 61 = 40° for asymmetric energy sharing of Ei : E2 ~ 
80% : 20%. Although not shown here, there is again 
good agreement between their results, Vanroose et al.'s 
one-center spherical ECS numbers [lll |. and our time- 
dependent FE-DVR predictions. 

It is also interesting to investigate the dominant es- 
cape modes for the various scenarios. These modes are 
strongly dependent on how the electrons share the excess 
energy. In an arbitrary geometry (0° ^ 0jv ^ 90°), for 
example, the back-to-back escape mode (6*12 = 180°) is 
forbidden for equal energy sharing. On the other hand, it 
becomes the dominant mode for significantly asymmetric 
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FIG. 7. (Color online) Coplanar TDCS of the aligned hydrogen molecule for asymmetric energy sharing. The electron detected 
at the fixed angle Oi — 90° takes away 20% of the available excess energy, while the second electron takes away 80% of i?exc. 
The present time-dependent FE-DVR results are compared with those from time-independent one-center spherical ECS [Tl|] 
and two-center prolate spheroidal ECS [11] calculations. 




FIG. 8. (Color online) Comparison of predicted relative coplanar TDCSs between H2 in the perpendicular (solid lines) and 
parallel (dashed lines) geometries, and He (chain lines) [42| for equal-energy sharing in polar coordinates. The polarization 
axis is taken along the horizontal direction. The photon energies for H2 and He are 75 eV and 99 eV, respectively. The fixed 
observation angles for one of the electrons are 0° (a), 30° (b), 60° (c), and 90° (d) with respect to the laser polarization vector. 
Scaling factors were used to emphasize the shape comparison. 



energy sharing, including the 20%:80% scenario discussed 
in the present paper (see Fig. [31). 

These results can be understood from a symmetry 
analysis [4l|. Equal-energy sharing and back-to-back 
emission is equivalent to fci = — fc2- When we con- 
sider the exchange and parity operations simultaneously 



in Eq. ((45|), we have ^-k^-ki = P(-l)'^$fci.fc2 • Here 
P = ±1 is the parity for the gerade and ungerade 
states, respectively. For the singlet double-continuum 
state with ungerade parity, we therefore must have 
^-k,k{'''i, = at any configuration of ri and r2- Al- 
though the magnitudes of the momenta k[ and k'2 are 
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FIG. 9. (Color online) Noncoplanar TDCS of the aligned hydro gen molecule for equal energy sharing. The present time- 
dependent FE-DVR results are compared with TDCC predictions [Ij]. 



not exactly conserved in the time-dependent picture, the 
ionization events we collect must satisfy the condition 
k[ = ^2 (because of the S function in Eq. ([55)) ) for the 
equal-energy sharing. This is the reason behind the for- 
bidden back-to-back {9i2 = 180°) escape mode for the 
equal energy sharing, as we observed in Figs. [5] and IB] 
Since the argument does not involve the relative align- 
ment angles, it is valid for all possible values of 9^. On 
the other hand, this is not the case when the excess en- 
ergy is not evenly distributed among the two electrons. 
Indeed, Figs.lSJa) and[3Jc) show maxima in the back-to- 
back emission, thereby illustrating the dramatic change 
in the dominant escape mode. 

For equal-energy sharing in the parallel geometry 
(6'jv = 0°), the electron-electron Coulomb repulsion sug- 
gests that the TDCS should be dynamically small if the 
two electrons travel along the same direction. This is 
in agreement with the numerically small cross sections 
(not exactly zero, however) at 62 = 0° or 360° seen in 
FigMd)- 

Recall that the one-photon double-photoionization 
process in helium [13, EH shares the same property. The 
back-to-back mode is forbidden for equal-energy sharing, 
and this can be explained by the above argument. It 
is one of the similarities between the molecular hydrogen 
and the atomic helium targets for double photoionization. 
However, Figs. [SI [51 and [71 also reveal significant molecu- 
lar effects in the TDCS results. These are missing for the 



helium atom, not only in the shape of the angular distri- 
butions, but also in the magnitudes of the cross sections. 
Depending on the relative orientation (0° < 6n < 90°), 
there is interference between the E„ and Hu symmetries 
in II2. A nice example of this effect was presented by 
Reddish et al. [13 . Even without interference (i.e., for 
6'jv = 0° or 90°), the perpendicular geometry shows much 
larger magnitudes of the TDCS than the parallel geom- 
etry. Figure [H shows the three cases of angular distribu- 
tions: II2 C -L e, II2 C II and He at equal energy shar- 
ing. Interestingly, in most cases the angular distributions 
of the perpendicular geometry resemble those of helium. 
The molecular effect can definitely not be ignored in the 
parallel geometry for 9i = 0° (c.f. Fig. [Slja)). The for- 
ward escape mode of the second electron is dominant for 
the H2 parallel geometry. In contrast, the backward mode 
is dominant for the H2 perpendicular geometry and also 
for helium. 

In Fig. [21 we show the TDCS for noncoplanar ge- 
ometries. Again, all angles are defined with respect to 
the polarization vector. For the perpendicular geometry, 
Fig. depicts the escape modes for the configuration 
of fci II C, (the fixed electron) and at the same time fc2 in 
the plane perpendicular to the plane formed by e and C- 
Figure shows the TDCS after exchanging the direc- 
tions of fei and k2 in Fig.lHlJa). With the same directions 
of fci and k2 as in Fig. [H^b), Fig. [5Jc) is for the case of 
the molecular axis orientated along the polarization vec- 
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tor. In the parallel case (^^r = 0°), we observe that any 
escape modes of both electrons ejected in the direction 
perpendicular to e arc forbidden. This can be under- 
stood by analyzing the spheroidal harmonics in Eq. (|53p . 
In this case, only the states can be populated. 

Hence, only partial waves with (mi, 7712) = {~m,m) 
and (—1)^1+^2 can contribute to the cross sections, since 
yeimiik', ki)ye^rn2{k' , ^2) at the angles of 6*1 = 6*2 = 90° 
vanish in spherical coordinates. Once again, the agree- 
ment between our FE-DVR noncoplanar TDCSs and the 
refined TDCC results [l3| is excellent. 



VII. SUMMARY 

We have presented calculations for one-photon dou- 
ble ionization of the hydrogen molecule at a photon en- 
ergy of 75 cV by solving the time-dependent Schrodinger 
equation in prolate spheroidal coordinates. The triple- 
differential cross sections were extracted through the pro- 
jection of the time-dependent wave packet onto uncorrc- 
lated two-electron continuum states, a few cycles of field- 
free time evolution after the laser pulse died off. 

Exhaustive convergence studies of the TDCS results 
were performed with respect to a number of discretization 
and expansion parameters, as well as the details of the 
laser field. These tests provide a strong indication that 
the results for the triple-differential cross sections pre- 
sented here are well converged and numerically accurate. 
Excellent agreement was obtained between the current 



time-dependent results in prolate spheroidal coordinates, 
those obtained with the ECS approach in spherical co- 
ordinates [m and, finally, larger TDCC calculations [l3| 
than those published earlier [13| . 

The present calculations do not confirm the signifi- 
cant reduction by about 20% in the TDCS results pre- 
dicted in recent ECS calculations in the two-center pro- 
late spheroidal coordinates [isj . Furthermore, our results 
did not show the level of sensitivity to the description of 
the ground state that was also reported by Tao et at. [l^ . 

The detailed analysis reported in this study provides a 
high level of confidence in the present results. We hope 
that they will be used as benchmarks for comparison in 
future investigations. Tables of these results are available 
in electronic format from the authors upon request. 
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